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We investigate the dynamics of Bloch electrons using a density operator method and connect this 
approach with previous theories based on wave packets. We study non-interacting systems with 
negligible disorder and strong spin-orbit interactions, which have been at the forefront of recent 
research on spin-related phenomena. We demonstrate that the requirement of gauge invariance 
results in a shift in the position at which the Wigner function of Bloch electrons is evaluated. The 
present formalism also yields the correction to the carrier velocity arising from the Berry phase. 
The gauge-dependent shift in carrier position and the Berry phase correction to the carrier velocity 
naturally appear in the charge and current density distributions. In the context of spin transport we 
show that the spin velocity may be defined in such a way as to enable spin dynamics to be treated 
on the same footing as charge dynamics. Aside from the gauge- dependent position shift we find 
additional, gauge-covariant multipole terms in the density distributions of spin, spin current and 
spin torque. 
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I. INTRODUCTION 

Carrier dynamics in metals and semiconductors in 
the presence of external electromagnetic fields, the po- 
tentials of which usually vary on scales considerably 
large than the interatomic spacing, have been conve- 
niently described by semiclassical transport theories. 
The semiclassical dynamics together with the Boltz- 
mann equation produce accurate descriptions of elec- 
trical and thermal conduct ionA. In the larger picture, 
semiclassical approaches are indispensable in problems 
involving both position and momentum, since in quan- 
tum mechanics position and momentum cannot be de- 
termined simultaneously. In recent years efforts have 
been made to extend the semiclassical theory to spin 
transport and generation 2 *^. The attempts to re- 
solve the challenges inherent in treating the transport 
of non-conserved quantities constitute a vibrant ongoing 
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A fundamental feature of semiclassical transport is its 
accounting for the finite extent of particles in real and 
reciprocal space. This feature is most naturally incorpo- 
rated into the dynamics of wave packet o 15 i 16 where the 
notion of a wave-packet center in real space and fc-space 
is retained. The carrier dynamics are described in terms 
of the displacement of these points under the action of 
external fields. The finite extent of the wave packet has 
important consequences for transport theory. For ex- 
ample, our recent research on spin transport has shown 
that, due to the fact that the spin and charge centers of 
a wave packet do not coincide, the expressions for the 
spin density, torque and current distributions are, in the 
language of wave packets, expressed as series of multipole 
terms 2 -. It has been demonstrated, in addition, that the 
wave packet formalism captures the physics connected 
with adiabatic motion and the Berry phase, in partic- 
ular the Berry-curvature correction to the semiclassical 
equations of motioni-V The Berry phase, which until 



Berry's seminal article^ 7 - was not taken into account, ap- 
pears naturally as part of the wave packet distribution 
function. The Berry-curvature correction to the wave 
packet equations of motion is believed to play an impor- 
tant role in the anomalous Hall eSec^^^ and in spin 
transports'* 4 **^, among other phenomena. 

Semiclassical transport theory is not restricted to wave 
packet dynamics. Wave packets, which may be con- 
structed out of one band eigenstatei 6 - or out of a super- 
position of eigenstates of several degenerate bands**** 2 *^, 
represent pure states. In order to treat mixed states (in- 
coherent superpositions of eigenstates) one must resort 
to a more general formalism. This is the principal moti- 
vation behind the current article. 

The most general description of a quantum mechanical 
system is based on the density operator. In this article 
we start from the density operator and formulate a the- 
ory of carrier dynamics in metals and semiconductors. 
We focus on non-interacting systems in which disorder 
is weak, strong spin-orbit interactions are present and a 
weak slowly-varying electric field is acting. These sys- 
tems have come under intense scrutiny in recent years 
along with the take-off of spintronica 23 i 24 i 25 i 26 i 27 . In such 
non-interacting systems the formalism may be simplified 
by defining a reduced one-particle density operatoi" 2 ^. 
Since the Bloch bands are clearly resolved the reduced 
density operator for this system may be expanded in a 
basis of Bloch wave functions. The density matrix which 
emerges from this expansion may be defined as a Wigner 
function. This function is used to study single particle 
dynamics and formulate definitions of macroscopic quan- 
tities. 

We focus on a number of fundamental aspects of adia- 
batic particle dynamics and demonstrate their relevance 
to transport phenomena. We pay particular attention to 
the band mixing induced by the electric field, which gives 
rise to a non-adiabatic correction to the wave functions. 
We also address several important questions concerning 
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the relationship between the Wigner function formalism 
and the wave-packet formalism^. We discuss the way the 
carrier position is to be found and, where possible, com- 
pare the result with the expression for the real-space cen- 
ter of a wave packet. The requirement that the particle 
position be gauge invariant results in a gauge-dependent 
shift in the position at which the Wigner function is 
evaluated. This shift was also found by Littlejohn and 
Flynn2£ in the study of coupled- wave equations. This 
gauge-dependent shift is a consequence of the freedom 
of choosing the real-space center of a wave packet by 
changing the phase of the wave functioni£*2i. The gauge- 
dependent shift in the position of the Wigner function 
must be taken into account when many-particle distribu- 
tions, such as the particle number density, are expressed 
in the crystal-momentum representation. The carrier ve- 
locity may be derived directly from the particle position, 
recovering the Berry-phase physics known from previous 
worki&. We also show that a spin velocity may be defined 
in such a way that spin transport may be described in an 
analogous fashion to charge transport. 

We discuss, in addition, important consequences of fi- 
nite particle size. The Wigner distribution is a quantum 
entity which takes into account the finite extent of the 
particles in real and reciprocal space. The distribution 
of, for example, spin for a single carrier, may not coincide 
with that of charge and the macroscopic spin distribution 
will be composed of a series of multipoles. The spin cur- 
rent and spin torque distributions are in turn composed 
of series of multipoles which we discuss and compare with 
those found in our previous work on spin transport 2 . 

The outline of this article is as follows. We present the 
fundamentals of the single particle density matrix formal- 
ism in section II. We determine the carrier position and 
velocity, emphasizing the gauge-dependent shift in the 
former and the Berry-phase correction in the latter. We 
calculate the particle spin, torque and spin velocity and 
define a modified spin velocity which satisfies an equation 
analogous to the charge velocity. In section III we demon- 
strate the modifications which must be made to extend 
the theory to many non-interacting particles. We define 
the charge and current densities and show the equation of 
continuity they satisfy. In the case of spin we define the 
spin, spin current and spin torque distributions and show 
the equation of continuity satisfied by them in the clean 
limit. In section IV we demonstrate the effect of local 
gauge transformations and the modifications which must 
be made to the dipoles in the charge- and spin-related 
distributions in order to make them gauge covariant. We 
conclude with a brief summary of our findings. 



II. SINGLE-PARTICLE DYNAMICS 

We consider systems described by a Hamiltonian hav- 
ing the following general form: 

H = H +H SO - (1) 



The term Hq is composed of the usual kinetic-energy con- 
tribution and the contribution due to the lattice-periodic 
potential. The term H so represents the spin-orbit inter- 
action term involving the carrier spins and the lattice- 
periodic potential. We restrict our attention to the limit 
in which this spin-orbit interaction is sizably stronger 
than the disorder broadening and the thermal broaden- 
ing. In this limit the system may also be described in 
terms of well-defined bands. The eigenstates of H , which 
have the periodicity of the crystal, are given by: 

H\ip m ) = e m \ip m ). (2) 

These wave functions are of the Bloch form, that is 
= e* k ' r |iii), where the functions \ui) represent the 
lattice-periodic parts of the \ipi)- The |ttj) are spinors 
with the full periodicity of the lattice. Since the Hamilto- 
nian contains strong spin-orbit interaction terms, which 
may depend on wave vector and position, it is not illu- 
minating to decompose the eigenfunctions into an orbital 
and a spin part. 

A. Density matrix 

We take the system under study to be described by 
a density operator p. It is not our concern in this 
article to determine an expression for the density op- 
erator for a given system, since methods of finding 
the density operator have been studied extensively in 
the past 32 i 33 i 34 i 35 i 36 i 37 i 38 i 39 , much work being done in 
the context of spin orientation of carriers in the non- 
degenerate limit. We assume the form of p to be known 
and study the role it plays in the dynamics of charge 
and spin. The density operator may be expanded in the 
Hilbert space spanned by a complete orthonormal set of 
Bloch wave functions \ip n k) as 

P = ^^/w( k ' k ')IV'™k}(VVk'|- ( 3 ) 

n,n' k,k' 

In the approach we follow, all the time dependence of the 
density operator is contained in the wave functions, so 
that p nn i(k, k') does not have time dependence. In ther- 
mal equilibrium it is diagonal and its elements, p„„(k, k), 
are equal to the Fermi-Dirac function fo(s n ), where e n is 
the band energy. 

The expectation value of any operator A is found from 
the formula 

(A) = TrpA = ^ ^ (^"k | P | Vv k' ) (ipn' k' | A | yj» k ) 

n,n' k,k' /a \ 

= tr^p(k,k')A(k',k). 

k,k' 

The notation A(k',k) stands for the matrix elements of 
the operator A, namely A n i n = (Vvk'^IVVik)- The op- 
eration denoted by tr is simply the matrix trace, which 
does not include the wave-vector summation. If the op- 
erator A is replaced by the identity matrix we obtain 
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Tr p = ^ k tr p(k, k). For a single particle, the normal- 
ization condition on the density operator is Trp = 1. 

In order to make transparent the analogy with the lan- 
guage of wave packets, center of mass and relative coor- 
dinates may be defined in fc-space such that Q = k — k' 
and q = k 1 2 k . The density operator as a function of 
these coordinates can be re-expressed as 



^2 ^2 Pnn ' ( q + ' q ) I ) ^ n 



(5) 



n.n' q,Q 



where \ip n ±) = e zq± ' r |it„-|-) and |u n ±) is the periodic part 
of the Bloch wave at q± = q ± ^ . 

The Wigner function corresponding to the one particle 
density matrix is found through the transformation 



/w(q+,q-) 



r e 



(q,r). 



(6) 



For the sake of concreteness the integrals are represented 
as three dimensional. Nevertheless, the theory applies 
to systems of any dimensionality. The Wigner function 
plays the role of a distribution function in the variables q 
and r. Technically, however, it is not a distribution since 
it may take negative values^. The inverse transforma- 
tion is 



Pnn' (q, r) 



Finally, replacing the vector summations by integrations 
we are able to represent the density operator in the fol- 
lowing form 



q+-(r-r) 



u n+ )(u n <-\e 



-iq--(f-r) 



(8) 



r 



where j dV — J J d ^yl ■ Fot the remainder of this sec- 
tion, the Wigner function p nn > (q, r) will frequently be 
abbreviated to p. The variables q and r are simply labels 
for the carriers, not physical observables. In particular, 
the dummy variable r in the Fourier expansion of the 
density matrix is simply the Fourier dual of Q and must 
not be confused with the position operator f appearing 
in the density operator. It does not correspond to an 
actual position. 



B. Carrier position 

If we consider a particle in one band, labelled n, p 
reduces to a scalar, which will be denoted by p n . The 
position of the particle can be found as the expectation 
value of the position operator, which yields 



Trpr= J dVp n {q,r)(r + K n ). 



(9) 



Here TZ n = TZ nn = {u n \i^^)- The integrand is not gauge 
invariant but the integral can be shown to be by chang- 
ing the variable of integration r to r' = r + TZ n . The 
connection TZ n has no position dependence, therefore the 
Jacobian of the transformation is unity. The expectation 
value of the position operator is 



Trpr = J dV p n (<i,r' -K n )r', 



(10) 



where dV' is defined in the same way as dV but with r 
replaced by r'. The gauge invariance of IjlQI) will emerge 



below. We conclude that r is to be interpreted as a label 
for the charge carrier, while the effective particle position 
is r' = r + TZ n . Neither the label r nor the gauge field 
1Z are by themselves gauge invariant, but together they 
form a gauge-invariant quantity which represents the true 
position of the carrier. This result was also found ear- 
lier, in somewhat different circumstances, in the work of 
Littlejohn and Flynn™. In the one-band limit, a clear 
connection can also be made with the dynamics of wave 
packets. The gauge-dependent shift in r reflects the free- 
dom of changing the phase of the wave functions \ipnk), 
I^Vi'k') in - It is t ne same freedom one has in defining 
the center of mass of a wave packet, demonstrated by 
Sundaram and Niui&, by changing the overall phase of 
the wave function^. 

For multiple bands, the expression for the particle po- 
sition is (the Einstein summation convention will be used 
henceforth) 

Trpf = J dV/w(q, r)(r£„/ n + TZ n > n ). (11) 
One can rewrite this expression as 



(12) 



Tr pi = tr / dV p(q, r) (r + TZ) . 



and make the substitution r' = r + TZ as in the single 
band case. Although TZ is a matrix the Jacobian of this 
transformation is unity. Finally, the expectation value of 
the position operator can be expressed formally as 



Tr/5f = tr / dV p(q, r' - TZ) r'. 



(13) 
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Unlike the single band case, the expression p(q, r' — 1Z) 
is a formal abbreviation for the Taylor expansion about 
r', that is p(q, r') — 1Z ■ V r /p(q, r') + higher order. 

Expression (|l(Jfl identifies the position of a particle. To 
determine if a particle is localized at its position one must 
calculate the variance of the position operator, the expec- 
tation value (f 2 ) — (f) 2 , and ensure that is does not di- 
verge. It is shown in the appendix that the variance does 
not diverge and that this result applies to any number of 
bands. 



C. Carrier velocity in an electric field 

We consider a system acted on by a weak external 
electric field. The effect of this electric field is incor- 
porated fully into the gauge invariant crystal momentum 
through the addition of the electromagnetic vector po- 
tential A = — Ei. Because the Bloch functions retain 
translational symmetry the electromagnetic vector po- 
tential does not enter the travelling- wave part of the wave 
functions, which have the form \ipnq) = e zq ' r |it nq ). The 
lattice-periodic functions |it nq ) depend implicitly on time 
only through the crystal wave vector q = q + ^ . How- 
ever, the presence of the external electric field results in 
a non-adiabatic mixing of the bands with the result that 
the perturbed lattice-periodic functions, |u mq ), have the 
following form to first order in the electric field^i 



l rnq) 



(14) 

The phase 4> m (q, t) includes the dynamical phase and the 
Berry phase. The differential is equivalent to q ■ A 

since J^|u„ q ) = 0. The result expressed by ifH)) is gen- 
eral. Moreover, although its derivation relies on the as- 
sumption that the bands are non-degenerate it can be 
shown that, when calculating intrinsic contributions to 
transport (for the definition of intrinsic please see the 
appendix and our recent work 3 ), the result also holds 
for degenerate bands with the difference that the sum 
must exclude all bands which are degenerate in energy 
with band m. The proof of this statement is given in the 
appendix. Henceforth, |u mq ) and \u m q} will be abbre- 
viated to \u m ) and \u m ) respectively. As stated above, 
given that the \u m ) are functions of q only, one may re- 
place ^ in by q • where q =— Equation IT^f) 
can then be written as 



\u m ) 



eE 



>«) , (15) 



where the connection TZ nm = (u n \i-^\u m ) . The \u m ) 
form a complete set. They are, however, not eigenstates 
of the time- dependent Hamiltonian H = H(q). 



In the evaluations of matrix elements in this paper the 
only property of the basis functions that is used is their 
Bloch periodicity. Therefore the results which are ex- 
pressed in terms of the lattice periodic Bloch functions 
hold as well for the \u n ) as for the \u n ). 

In the absence of disorder, since all the time depen- 
dence is contained in the wave functions, the density 
matrix j o„„'(k, k') in @ can only depend on the wave 
vector k, not on the crystal wave vector k. As a result 
the Wigner function p nn i (q, r) only depends on the wave 
vector q, not on q. 

It is customary to consider only a subset of the Hilbert 
space which contains the bands that are relevant to trans- 
port, which in semiconductors usually refers to the top- 
most filled valence bands and/or the lowest filled conduc- 
tion bands. Since the gauge of the \u m ) is not fixed we 
impose the following gauge-fixing condition in the sub- 
space under consideration 



u n \ih-^\u m ) — (u n \Ti.\u m }, for n in subspace 
0, for n out of subspace, 



(16) 



where r fi = e~* qr _ffe lqr . This condition fixes the 
phase(s) 4> m in l|14|) . We shall henceforth work only with 
the basis set {|u„)}. 

The particle velocity can be derived directly from 
by evaluating the time derivative 



d m " f n, dH 

-Tr P r = tr dV P -. 



(17) 



Using equation (|l(j[) . the differential becomes, after a 
little straightforward algebra, 



1 dH n > n . / du n > du n du n > du n , 



(18) 



The abbreviation H n > n stands for the matrix elements of 
the time-dependent Hamiltonian H in the presence of an 
electric field, (u n >\H\u n ). In one band it is easily shown 
that (|T%)) becomes 



dt ~hdq q X S *" + -« ■ 



(19) 



The Berry curvature for band n, fi n , is given by 
V q x TZ n and V q x represents the curl operator in 
reciprocal space. The quantity H^ q is defined as 

* " <n^W)- Equation © is the ana- 

log of the semiclassical equation of motion found by Sun- 
daram and Niui£. Writing q =— Tp we obtain for elec- 
trons in a single band: 



(20) 



For multiple bands, an elegant result is obtained by in 

<)q 



troducing the covariant derivatives^ = — i1Z and 
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Charge 



Spin 




FIG. 1: For a particle of finite extent the charge and spin 
distributions in real space are in general do not coincide. The 
same is true of the charge and spin distributions in reciprocal 
space. 



_D 

Dt 



d_ 

<// 



iA, where A n >n = ( u n'K^ L )- Making use 
of these derivatives, equation (|18|l can be written in the 
manifestly gauge covariant form 



dlZ a , I 



dt 



D 



h Dq a 



D D 



L Dq a ' Dt 1 



(21) 



It is important to point out that, from the gauge- fixing 
condition which defines A, it is evident that A is 
gauge covariant, implying that the time derivative ^ is 
itself gauge-covariant. This is a peculiarity of the gauge- 
fixing condition we have chosen. 



FIG. 2: In the presence of spin-orbit interactions the spin 
distribution of a particle changes in time. The horizontal axis 
may represent position or wave vector. 



distributions may be quite different. These facts are il- 
lustrated in Figs. 1 and 2. 

To evaluate the average carrier spin consider the ex- 
pectation value of the operator s, which stands for any 
one component of the spin operator. The result is 



Tr ps = tr / dV ps, 



(22) 



where s„'„ = (u n i\s\u n ) are the matrix elements of the 
spin operator. In the presence of spin-orbit interactions 
a spin torque is associated with the spin of every carrier, 
which accounts for the non-conservation of spin, as shown 
in Fig. 2. The average carrier torque is found by evalu- 
ating the expectation value of f = i[H, §]. The result is 



D. Carrier spin, spin torque and consequences of 
finite particle size 

The fact that carriers have a finite extent in real and 
wave-vector space has profound implications for parti- 
cle dynamics and transport. These implications were 
pointed out in a previous publication^ using the semi- 
classical language of wave-packets and will be elaborated 
in this section from the density matrix point of view. For 
the sake of clarity we will specialize in spin, although the 
discussion in this section applies to any quantity. 

When considering the transport of spin in a system 
composed of many carriers, one must associate a spin 
distribution with each individual carrier. A center of 
spin may be defined for each particle by r s = This 
definition may be troublesome if the expectation value of 
the spin operator were zero, but its only purpose is to 
illustrate a physical principle. Evidently, if one replaced 
the spin operator with the charge or mass one would ob- 
tain the particle position as given by ©, which we will 
refer to as r c . It is obvious from the definition of r s that 
in the general case there is no reason for the center of 
the spin distribution of an individual particle to be the 
same as the actual particle position. This center will be 
different for each component of the spin operator. Fur- 
thermore, since spin is not conserved in the presence of 
e.g. spin-orbit interactions, r s may also be a function of 
time. This suggests that the spin distribution of one car- 
rier will in general have a different shape than its charge 
distribution, and that the time development of the two 



Trpf = tr J dVpf, (23) 

where f n > n = (m„'|t|u„). 

As a result of the distinction between r s and r c , in cal- 
culations of spin-related quantities which use the center 
of charge as the reference, multipole terms must be taken 
into account in addition to the average quantities calcu- 
lated above. For example, in the distribution of spin a 
spin dipole will be present, as well as higher order mul- 
tipole terms, which are assumed small. The spin dipole 
is found as the expectation value (rs), evaluated in the 
appendix, and yields the gauge field 

M^ n = i((^|J|^)-(^|JK)). (24) 

In contrast to r s , M s is well defined. Similarly, in the 
distribution of torque of an individual particle a torque 
dipole will be present. The torque gauge field M r = (ff) 
is given, as shown in the appendix, by 

Ml, n = i((u n ,\r\^)-(^\r\n n )). (25) 

M s and M T are gauge dependent. It will prove useful 
to define gauge-covariant dipoles as p s = M s — ^{1Z,s} 
and p r = M r — ^{TZ, ?} respectively. The remainder of 
the discussion of these gauge-covariant spin and torque 
dipoles is deferred to section IV. 
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E. Carrier spin velocity in an electric field 

Wc define the spin velocity as the expectation value 
(sv), where products of non-commuting operators are 
assumed to be symmetrized. The result is 



Tr psv = tr / dV p v 



(26) 



The spin velocity v s is given by 

d 

v s = — 
dt 

dM 



(f s) — (rf ) 
-M r . 



(27) 



dt 



In order to write the velocity in a gauge covariant form, 
we replace in the above M s = p s + \{R, s} and M T = 
p T + 7j{it,?}. This yields immediately 



i <m dp s 
2 { ^' s| + ^T 



(28) 



Equation l|28|) is the gauge-covariant form of equation 
(10) in Ref. [2], the integrand of which represents the 
spin velocity. The first term in (|28|) is a convective con- 
tribution and represents a moving electron transport- 
ing its average spin along with it. The second term 
is the time derivative of the spin dipole while the last 
term is the torque dipole which takes into account the 
non-conservation of spin. We may incorporate p r into 
a modified spin velocity 2 *^, which we shall call v*, by 
v = v s + p T . The modified spin velocity is given simply 
by: 



1 dK 
2^~dt' 



s} + 



dp s 
~dt 



(29) 



The importance of this definition of the velocity will be 
seen in the following section, when the macroscopic spin 
current is introduced. 



III. MANY PARTICLE DISTRIBUTIONS 

We will consider the macroscopic distributions of 
charge and spin starting from the formalism we have de- 
veloped. In particular, our discussion of the corrections 
to the spin density, current and torque distributions is 
motivated by the observation that the existing literature 
has omitted various contributions to spin transport. Sev- 
eral works have used semiclassical concepts but did not 
arrive at answers containing all the terms we have de- 
rived. 

The particle number density is defined by the formula: 



n(R, t) = Tr[p<5(R - f ) 



(30) 



The trace operation Tr involves integration over r. The 
procedure we follow has affinities with the coarse graining 



of electrodynamics in material medial. The size of the 
carriers is taken to be smaller than the length scale of the 
Wigner function. We regard the 5-function as a sampling 
function^ with a width somewhere between the micro- 
scopic scale of the carriers and the macroscopic scale of 
the Wigner function. The rationale for regarding the 
(^-function as a sampling function comes from the real- 
ization that often the physics of the problem does not 
require absolute resolution of position, provided that the 
resolution is finer than the length scale of the external 
field. Since the variance of f is finite, the (5-function may 
be expanded in a Taylor series about the dummy variable 
r in the form 5(K- f ) = <J(R-r) - Vr • (f - r)<5(R - r) + 
0[(r — r) 2 ]. The expansion is truncated at the first order 
for simplicity but we will show in the next section that 
all the terms may be recovered in a concise and elegant 
fashion. The number density can be expressed as: 



n(R,t)=tr / d 3 q(p-V R - pK). 



(31) 



In the above p stands for p(q, R), an abbreviation which 
will be frequently used in the remainder of the article. 
Note that the gauge field TZ in (jBTjl plays the role of a 
dipole. 



A. Electrical charge and current densities 

The charge density is defined by the formula 

n q (R,t) = qn(R,t) = 9 Tr[p<5(R - f)]. (32) 

The charge q is not to be confused with the wave vector q. 
The charge current density J 9 is defined by the equation 



J 9 (R,<) = qTr[pS(R - r)ir}. 



(33) 



The charge equation of continuity in the absence of ex- 
ternal sources and disorder, 



<9n« 
~dt 



V • J 9 



0, 



(34) 



is readily verified from the first-principles definitions of 
the charge and current densities. 

When the (5-function is expanded the current density 
takes the form 



J 9 (R,t) = qti 



d 3 q 
(2^)3 



•pv. 



(35) 



The velocity matrix elements v„<„ = {u n i |v|m„) as shown 
in the appendix, where v = e~ tq ' r ve lq ' r . Since in the 
equation of continuity it is the gradient of the current 
that appears, we have not included in corrections to 
the current which arise from the fact that the velocity 
distribution of a single carrier is different from its charge 
distribution. These corrections are in principle present 
but have been omitted for simplicity. 
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In homogeneous systems the gradient term in (|35|) 
drops out of the equation of continuity and, based on 
equations Ijl7|l and the charge current can be ex- 

pressed as the time rate of change of the gauge field TZ 

B. Spin, spin current and spin torque densities 

The spin density is defined as 

5(R, t) = Tr[/55(R - r)s], (37) 

while the torque density takes the form 

T(R, t) = Tr[p5(R - r)f] (38) 

and the spin current density is defined as 

J s (R,t) = Tr[p<5(R-r)sv]. (39) 

In the absence of disorder the spin equation of continuity 
is 

— + V • J s = T, (40) 

which is verified from the first-principles definitions in- 
troduced above. 

The (^-functions are expanded in the same way as for 
the particle number density, whereupon the spin density 
can be expressed as 

S(R, t) = tr J -0^ (p S "-Va-pM s ), (41) 

and the torque density is 

T(R,i)=tr|^(pf-V R .pM T ). (42) 

Ignoring the gradient term, the spin current is 

JS ( R '*)= tr /(0^ S - ( 43 ) 

In analogy with the modified spin velocity of the previous 
section, a modified spin current J* may be defined by: 

J ^ (44) 
f d 3 q (\ c dK dp s \ v ' 

The equation of continuity satisfied by this current in the 
clean limit is: 

ds „ Tt f d3 <i 

m +v - J - tr 7 Wf pf - (45) 



In many models, such as the spherical four-band Lut- 
tinger model^, the RHS vanishes when all bands in the 
subspace are taken into account. In that case the spin 
current J* is conserved. Much effort has been devoted to 
finding a conserved spin current. In our previous work»2ii£ 
we have argued, based on semiclassical ideas, that the 
closest one can come to a conserved spin current is by 
including the torque dipole in the definition of the cur- 
rent. In this article we have derived this current from a 
density matrix point of view and show that it supports 
our earlier conclusions. To date, the definition J* is the 
closest one has come to a conserved spin current. More- 
over, as shown also by Zarea and Ulloa™, the behavior 
of the modified spin current may be rather different from 
that of J s . 



IV. GAUGE TRANSFORMATIONS 

We have argued in section II that in discussing trans- 
port of an observable one must take into account the fact 
that the distribution of that observable for each individ- 
ual carrier in general involves a dipole correction, and in 
principle also higher-order corrections. If the quantity 
being transported is not conserved then the determina- 
tion of its rate of change also involves a dipole correction. 
The existence of what appears to be a sound physical 
argument for the inclusion of these corrections suggests 
that objects such as the spin dipole and the torque dipole 
are not simply mathematical artifacts meant to facilitate 
calculations. It is appropriate to investigate whether they 
are in fact fundamental quantities with a true physical 
meaning. If that were the case, they ought to be ex- 
pressible in a gauge-covariant way and to give rise to 
observable effects. 

In addition to these considerations, our aim of provid- 
ing a physically transparent formalism requires a test of 
the gauge covariance of the macroscopic quantities and 
equations of motion we have derived in this article. The 
physics contained in them must be independent of the 
choice of basis functions. Furthermore, the physical dis- 
cussion can become cluttered with futile and meaningless 
objects if it is formulated in terms of gauge-dependent 
objects. In contrast, the gauge-covariant expressions we 
derive below are simple, elegant and transparent. 

A general local gauge transformation is represented by 
\u m ) — > O mn \u n ), with O mn = mn (q) and m and n 
are indices of bands in the subspace under consideration. 
Under this operation 

BO 

K^K + iiO- 1 —), (46) 
aq 

in which TZ stands for OTZO^ 1 . The Berry curvature for 
one band O n is discussed extensively in the paper of Sun- 
daram and Nivj^. One important additional detail which 
emerges from the transformation given by (|46|l and the 
definition fl n = Vq x 1Z„ is the fact that, if the curvature 
is non-zero, one cannot make a gauge transformation to 
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eliminate K n . It can be easily shown that, whereas K n is 
gauge dependent, the Berry curvature is gauge covariant. 
Therefore if f2„ is non-zero in one gauge it is non-zero 
in all gauges and there is no gauge in which K n can be 
zero. As a result, in systems in which the curvature is 
non-zero the gauge-dependent position shift is necessar- 
ily present in the Wigner function. It can be regarded 
as a 'penalty' for working with the position operator in 
a basis of definite wave vector. 



n(r) 




Real position 



Apparent position 



A. Gauge covariance of macroscopic densities 



FIG. 3: The position which appears in the Wigner function 
is not the real position. 



The Wigner function itself changes under the gauge 
transformation defined above. Using ©, we find that, 
to first order in O, the Wigner function changes as 



(47) 



where p = O 1 pO (the opposite of the gauge connection 
K). 

All the macroscopic densities defined in the previous 
sections are covariant under a local gauge transforma- 
tion. This will be demonstrated for the spin density. 
Under the above gauge transformation, the gauge field 
M s transforms as 



M s 



(48) 



The matrix elements of the transformed gauge field M s 
are given by M s n , n = § ( (u n ,\S | % ) - ( %i 1 5 \ u n ) ) and the 
transformed spin operator s = _1 sO. To first order the 
extra terms acquired by the spin density under a gauge 
transformation are 



BO 

{^O- 1 — }i)=0, (49) 



so that the spin density remains gauge covariant. The 
cancellation remains true for all orders in the expansion. 
Similarly, the charge and current densities do not acquire 
additional terms under the local gauge transformation in- 
troduced above. The extra terms appearing as a result of 
the transformation cancel when the trace is taken. When 
the change in the Wigner function under a gauge trans- 
formation is taken into account the overall expressions 
are gauge covariant. 



B. Gauge-covariant expressions for spin and torque 
dipoles 

Equation l|31[) for the particle number density can be 
formally written the following way 



rc(R,t)=tr / d^qpi^R-K). 



(50) 



Re-expressing the integrand in this manner is tanta- 
mount to making, in the density matrix, the replacement 
p(q, r) — > /?(q, r — K). To ensure the gauge covariance 
of the number density the Fourier dual r is replaced by 
the true position r + K, as illustrated in Fig. 3. This is 
the same position as found in section II. If the subspace 
contains one band expression (|50(l is an exact result, not 
simply a formal way of writing the number density. 

Examining more closely the spin density as given in 
(|41|l it is evident that, whereas the spin density it- 
self is gauge covariant, its individual constituents are 
not. If one were to consider the integrand of the dipole 
term in (|41|l without the Wigner function, this quan- 
tity would not by itself be gauge covariant. We have 
already constructed a gauge-covariant spin dipole p s = 



-,{K, s}. In terms of the gauge-covariant spin 



dipole, the spin density can be re-expressed as 



S(R,t) =tr 



d 3 g 
(2tt) 



(ps- V R -pKs- V R -pp s ). 



(51) 



It can be formally written the following way 



S(R,t) = tr 



rf 3 



q 



[p(q, R-K) s-V R -p(q, R-K) p']. 

(52) 

The gauge-covariant spin dipole p s , as pointed out above, 
is the result of a carrier's center of spin being different 
from its center of charge. 

In the same way the gauge-covariant torque dipole has 
been defined by p T = M T — \{TZ, f}. Carrying out an 
identical manipulation to that for the spin density, use 
of the gauge-covariant torque dipole allows us to rewrite 
the torque density formally as 



T(R, t) = tr 



d 3 q 
(2tt) 3 



[p(q, R-K) f-V R -p(q, R-K) p T ]. 

(53) 

A similar formal expression exists for the spin current. 
Evidently, since these expressions are simply a rewriting 
of the spin, torque and current densities, the equation of 
continuity satisfied by them is unaltered. 

We remark that the gauge-covariant quantities intro- 
duced throughout this paper are well defined. In fact, 
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even though the density matrix was expressed in a basis 
of extended Bloch states, by means of various manipula- 
tions we have been able to obtain well-dehned formulas 
for all the objects relevant in transport. 

Finally, as described by Dudarev et alm& in the con- 
text of atomic physics, optical lattices can be constructed 
which mimic the spin-orbit interaction. In such a lattice, 
a wave packet can be constructed using cold atoms and 
the evolution of its centers of mass and spin can be fol- 
lowed. In the same spirit, Kato et a/4£ have shown that 
it is possible to follow the motion of a spin packet in 
solid state systems, which demonstrates the feasibility of 
an analogous experiment in InGaAs. This would provide 
a way to measure a spin dipole directly. 



For a set of degenerate bands the equilibrium part of the 
density matrix is proportional to the identity matrix / 
and will have the form /o /, where /o is usually the Fermi- 
Dirac distribution. Transport theory often distinguishes 
between intrinsic effects, which are due to the equilibrium 
part of the density matrix, and extrinsic effects, which are 
due to the non-equilibrium correction to it. If we wish 
to evaluate the expectation value of an operator A using 
the equilibrium part of the density matrix for a set of 
degenerate bands, which for simplicity here we will take 
as being two dimensional, the following quantity must be 
evaluated: 

(A) = fodiLilAlu!) + (U 2 \A\U 2 )) 



V. SUMMARY 

We have formulated a theory suitable for describing 
the dynamics of particles in both single and multiple 
bands. The one-band results known from the wave- 
packet formalism, including the terms connected to the 
Berry phase, emerge from our theory. The formalism 
can be applied to any clean system regardless of the di- 
mensionality of the Hilbert space under consideration. 
As a result of the gauge degree of freedom of the basis 
functions, the position vector which is used as a label of 
the carrier must be modified by a gauge-dependent shift 
in order to obtain the true particle position. We have 
shown the way to define macroscopic density distribu- 
tions for conserved and non-conserved operators. In the 
case of spin we have highlighted the correspondence be- 
tween the multipole terms which appear in the density- 
matrix formalism and those found earlier in the wave- 
packet formalism and we have shown that experiments 
can be performed to identify the effect of a spin dipole. 

DC was supported by the NSF under grant number 
DMR-0404252. QN was supported by the DOE under 
grant number DE-FG03-02ER45958. 



VI. APPENDIX 



The perturbed wave functions are given by: 



\ui) = \ui) - 



\u 2 ) = K) 



(U2\ih4f\ui) ^ (UnliHjj-Aux) 



Si - £2 

( Ul \ihf t \u 2 ) 

£ 2 - £l 



"2 



- m £ n 



«1 



{u n \ihj- t \u m ) 



u„) 



\u n ). 



n/1,2 



The expectation values are: 



(ui\A\ui) = A n - 



; du-2 
d±_ 

£\ - £2 

dU2 

dt 



A Hjffi 

^21 



£\ - £2 



■A 12 + J2 



1 

out 



(u 2 A\u 2 ) = A 22 - { - 11 dt 1 A 21 - V 21 dt ' A 12 + V 
£2 ~£\ £2 - El 



In the above Y]° u and ^2" stand for the sums involving 
the bands outside the degenerate manifold and Ay = 
(«,|A|uj). The denominators of the terms involving A\ 2 
and A 2 \ have opposite signs. Therefore, adding up the 
expectation values we obtain: 

out out 

«ui|A|tii) + (u 2 \A\u 2 )) = A n + A 22 + ^ + . 



We will present in this appendix some of the proofs 
and evaluations which require lengthier calculations and 
would interrupt the flow if incorporated into the main 
text. For all derivations except the first, the \u m ) and 
\u m ) are interchangeable. 



Therefore the terms with diverging denominators cancel 
out. 



B. Evaluation of position matrix elements 



A. Proof of Eq. (1141) for degenerate bands 

The perturbed eigenfunctions are given by Ea. (|14f> : 

v (u n \mj\n m ) 



In the general case, the expectation value of the posi- 
tion operator is given by 



f) = / dV p nn > d A Qe 



d 



dV Pnn , I d 3 Qe-^ r (u n ^\(^i—e l ^)\u n+ ). 
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We can expand the integrand in the above equation and 
obtain 



{u n ,-\ie t(: *- i \u n+ ) = -i — (u„/_|e iQ ' f \u n+ ) 



,du n i- 



e iQ 'W> 



+ <^-K Q - f |^> 
+r(u„/_|e iQ ' r |u„ + ). 



q, treating Q as a small parameter 

Q i du n£ . 



\u n+ ) = |u„q) + — ' 



(54) 



2 '^q" 7 

|«n-> = l«nq> - y ' l^p) 

In the limit in which Q — > 

1 9Q ' ~ 2 1 9q 

Consequently 

(u n '— re ^ " u n +) rt5 n / n 7^. n / n , 

and the expectation value of the position vector yields, 
finally 

(r) = y dVp nn /(r5„/ n + H n >n)- 

C. Evaluation of spin and torque gauge fields 

The spin gauge field is found in an analogous fashion 
to the expectation value of the position operator 

(rs)= J dV Pnn . J d 3 Qe-^- r (S„^|fse lQ - f |S rl+ ) 

d_ 
d~Q' 

= [ dV Pnn ,M s n , n . 



= J dV Pnn , J d ■-<■>: '■«■*•{»„, |.s(-/Ji : r'-Q-'-)| ( ,„_ 



Expanding the above, 



d 



(u„/_|fse 4Q - f |u n+ ) = -i-^(u n ,-\se lQ ' r \u n+ ) 
+{ u n ,_\se^\i^). 



Using the same arguments as in the previous section, 
the first term integrates to zero, and after evaluating 
the brackets involving the exponential and the lattice- 
periodic functions, we obtain in the limit Q — > 

W> n = % tiM*\^)-(^\*K))- (55) 

An almost identical derivation applies for the expecta- 
tion value (rf) = ^({f,f}}. The only difference is that, 
since f may be a function of wave vector, an additional 
term involving ~i§^ is generated. However, that term 
drops out when the anti-commutator is taken, leaving us 
with the result 



M n , n = -((u n ,\T\—)-(— 



\r\un))- (56) 



D. Finite variance of carrier position 

Since the expectation value of the carrier position has 
been shown above to contain no divergences, if the ex- 
pectation value (f 2 ) is not divergent then the variance 
(f 2 ) — (f) 2 is also finite. To see that the variance of 
the particle position does not contain any divergences, 
we need to evaluate the expectation value (f 2 ). This re- 
quires us to evaluate the matrix element 



(v 



U n+ ) = -(tt„/_|( 



dQ 2 



)\Un+). 



First note that an expression of the form —A^^CD can 
be written as 



-A^CD 



-w [ABCD) 



d 2 A 



BCD 



dQ dQ K ' 



In the case we are considering, all the products involve 
brackets which are proportional to 8{Q). As a result, all 
the terms on the first line vanish under integration with 
respect to Q and only the terms on the second line need 
to be evaluated. Writing them out explicitly, 



(u n ' 



|fV c 



\u n+ ) = r 2 (u„,_|e iQ - f |u„ + )) 



dQ 2 
+2ir-(( 



\e^\u n+ ) - (u n ,.\e iC >- f \ 



d 2 u n+ v 
dQ 2 1 



du n >- 
dQ 



e^\u n+ ) + {u n ,_\e* 



dUn^ 

1 <9Q 



» 



-2 



dQ 



du Tl 



dQ 
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In order to evaluate the differentials, the wave func- 
tions are expanded as in (|54|l except that now the expan- 
sion must be made to second order in Q. The final result 
is 

(f 2 ) = J dV p nn <[r 2 6 n . n + 2r • TZ n , n 
1 du n , du r 1 d 2 u n > 

d 2 u r . 
dq 2 



This is clearly finite so the variance of the position oper- 
ator is finite. 



All the brackets in the above two equations are propor- 
tional to <5(Q), causing most of the terms to cancel. We 
are left with: 



1 



V-e 



— iq_ ■ (r— r) 



{v,(f-r)}e«i+^- r )K + ) 



^(Q)««»m^>-<^ivi^». 



F. Effect of gauge transformation on p(q, r) 

The operator p must be invariant under gauge trans- 
formations. From the expansion of the density operator 
in Bloch eigenstates, 



E. Evaluation of velocity matrix elements 

The velocity matrix element (u n >— |e~ ,q ~' r ve zq+ ' r |£t n +) 
is easily evaluated as 

= J d 3 r' v_(r>- , i- r ve !q+ ' r 's„ + (r') 

= J d 3 r' e jQ - r 'w n /_(r')(e- iq +- f 've jq +- f ')u„ + (r') 

= 5(Q) J d 3 r'u n ,(r')vu n (r') = 5(Q)(u n ,\v\u n }. 

In the above we have abbreviated v = e~ iq ' r ve lq ' r . 

We also need to evaluate the commutator 
i(u„,_|e- 4q --( f - r ){v, (r - r)}e iq +-^ r )|u„+). One 
half of the commutator is evaluated as 



-iq— ■ (r — 



v(f — r)e 



iq + -(f— r) 



U n+ ) 



-iq_-(f-r) • 



. d 

<9q 4 



D iq + -(f-r)i 



Wn+) 







= -Itt- [ u m - 1 e" 4q - ■ (f ~ r) ve iq + u n+ ) 
dq+ 

oq + 

+ (^_|e-' iq --( i - r )ve 4q +-( f - r )|i^±), 



dq+ 



while the other is 



V- e 



iq--(r-r)/- _ r j^ e «q+-(f-r) | 



(u n ,_\\i^e- t < 1 -^- r) }ve^< i -^\u n+ ) 



<9q_ 
,. du n - 
9q_ 



dq- 

| e -iq_-(f-r)« iq+.( 



VI - r ^ r> \u n+ ) 



iq--(f-r) - iq+-(f-r 



' " ' ■v<-' ■ ■'\u„> + ) 



{u n '- e 



iq--(r-r) OW c »q + -(r-r) 

9q- 



n,n' k,k' 

it is evident that, if the wave functions change according 
to l^mic) OmnO^l^nk), then yO n „'(k, k') must trans- 
form to 0~ 1 pO. Using the definition of p nn >(k, k') in 
terms of the Wigner function, 

P(q,r)=y ^-^e jQr p(q + ,q_), 

and remembering that k = q + and k' = q , we obtain 
d 3 Q 



p( 



q,r) - J ^ e i Q- r O- 1 (q + )p(q + ,q_)0(q_). 



The matrices O 1 (q+) and O(q-) may be expanded 
about their arguments as: 



0- 1 (q+) = 0- 1 (q) + 



Q dO-^q) 
2 ' dq 



0{Q 2 ) 



0(q-)=0(q)-§- ^f + 0(Q 2 ) 



Then, to first order in Q, abbreviating 0(q) to O and 

p(q+,q-) to p g , 

0- 1 (q + )p g (q +) q_)0(q_) 



where p q — O 1 p q O. The last line is obtained by insert- 

Jq^ 



ing OO 1 or O 1 as appropriate. Since -§^{00 1 ) = 0, 



we have ^-O = -Q- 1 ^- and 



1 dO 



9q 



. Q s ~ n -idO 



Finally writing — Q = |V r e i( ^' r , we recover formula l|47|) 
for the transformation of p(q, r). 
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